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We study the natural excitation of molecular systems, applicable to, for example, photosynthetic 
light-harvesting complexes, by natural incoherent light. In contrast with the conventional classical 
models, we show that the light need not have random character to properly represent the resultant 
linear excitation. Rather, thermal excitation can be interpreted as a collection of individual events 
resulting from the system’s interaction with individual, deterministic pulsed realizations that con¬ 
stitute the field. The derived expressions for the individual field realizations and excitation events 
allow for a wave function formalism, and therefore constitute a useful calculational tool to study 
dynamics following thermal-light excitation. Further, they provide a route to the experimental 
determination of natural incoherent excitation using pulsed laser techniques. 


“One of the most important attributes of a 
stationary random process is its spectrum.” 

L. Mandel & E. Wolf [1] 

The process of light-matter interaction can initiate fun¬ 
damental dynamics in systems as diverse as semiconduc¬ 
tors, quantum dots, and living organisms. Of partic¬ 
ular interest is the linear interaction of sunlight with 
molecules such as photosynthetic light-harvesting com¬ 
plexes [ 2 ] , which is the foremost step responsible for pri¬ 
mary energy production on earth. 

While the development of advanced experiments 
such as non-linear spectroscopy [3] allows characteriza¬ 
tion of the system’s dynamical properties upon photo¬ 
excitation, the measured molecular dynamics can differ 
from that resulting from natural excitation conditions 
[4-9], mainly because of significant differences between 
the exciting lights [10]. Although in theory, it is possi¬ 
ble to measure molecular dynamics directly from inco¬ 
herent excitation [ 11 ], technical difficulties are such that, 
in practice, no experimental data are yet available. At 
present, the dynamics of a system excited by thermal 
light can only be inferred from appropriate theoretical 
models, which require a proper description of thermal 
excitation. Two different approaches are available: (i) 
a full quantum-mechanical treatment, using a quantum 
description of the light (well-known for thermal light, e.g. 
the P-representation, representing the light state by di¬ 
rect products of incoherent mixtures of coherent states 
associated with a single frequency—typically plane waves 
[ 12 ]), or (ii) a semi-classical approach—based on a clas¬ 
sical description of the radiation field. 

Thermal light can only be represented statistically by 
a probability distribution function or, alternatively, as a 
set of all possible realizations. The most conventional 
classical representation - collision-broadening model [ 12 ] 
- relies on realizations with stochastic phase jumps to 
reproduce the random character of thermal light, with 
the average time interval between jumps reproducing the 
(short) coherence time. Here we show that no random 


parameter is needed for the proper representation of lin¬ 
ear excitation, i.e. weak excitation, and that, in first- 
order perturbation theory, an ensemble of deterministic 
realizations can reproduce the result of thermal-light ex¬ 
citation. Indeed, even pulses are possible. In studying 
the interaction of a molecule with an ensemble of helds 
(representing the individual realizations), we provide an¬ 
swers to the following question: What field characteris¬ 
tics provide a proper description such that the ensemble 
of individual realizations yield the correct molecular den¬ 
sity matrix? 

We address this question within the framework of a 
molecule, applicable to, e.g., photosynthetic pigment- 
protein complexes such as the Fenna-Matthew-Olson 
(FMO) complex [13], or the LH2 light-harvesting com¬ 
plex [14]. Because we focus on the molecular response to 
the light, the relevance and interpretation of our results 
are independent of the effects of the protein environment, 
and we assume a closed molecular system. 

Consider the Hamiltonian 

H{t) = Y,e^\a){a\-d - {t) + (t)*) , ( 1 ) 

OL 

where the system Hamiltonian (first term on the r.h.s.) 
is written in the basis of eigenstates [{a:}) with eigenen- 
ergies and the second term describes the interac¬ 
tion with the field in the semi-classical approach. Here 
d = Ta\oi) (ffjnQ+b.c., with \g) the molecular ground 
state, is the transition dipole moment operator, oriented 
along the unit vector n^; E(''')(t) denotes the positive 
frequency part of one individual field realization. The 
superscript 7 is used throughout the paper to label a 
particular realization, either of the held, or of the cor¬ 
responding molecular excited state. This Hamiltonian 
characterizes the molecular dynamics upon interaction 
with one of the held realizations. 

In the rotating wave approximation and in the limit 
of weak-held interaction, hrst-order perturbation theory 
provides a solution for the wave function resulting from 
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the \g) state excitation as a superposition of eigenstates: 

a=l do 

( 2 ) 

where Na is the number of molecular eigenstates and 
uia denotes the eigenfrequencies. Here, we assume the 
excited-state level splitting to be larger than their spon¬ 
taneous decay widths and neglect effects from radiative 
decay [ 15 ]. For the sake of clarity, we consider a one¬ 
dimensional field, and assume it lies along the orientation 
of the transition dipole moment ha, taken identical for 
all eigenstates. 

Considering that the molecule’s response to a stochas¬ 
tic field is composed by an ensemble of individual exci¬ 
tation events, we define the density matrix: 

( 3 ) 

where (.. .)p denotes averaging over the ensemble of re¬ 
alizations. From ( 2 ), the elements of this density matrix 
are: 

= J (ir2 (4) 


with OJaP = OJa — OJp. 

Now we require that this density matrix be represen¬ 
tative of thermal-light excitation, and obtain the restric¬ 
tions on the fields and ensemble (.. .)p. For the 

molecule considered in (1), interaction with thermal light 
would lead to the excited state described by the density 
matrix with elements: 

f f dr2 

^ Jo Jo ( 5 ) 

X -T2). 


To calculate the ensemble-averaged auto-correlation 
function of the field, we use the Fourier representation 
for the individual realizations [22]: 

pOO 

E('r){t)= duj e‘'^\uj) ( 8 ) 

Jo 

where E^'^\u)) denotes the Fourier components. We 
choose all individual fields with identical spectral distri¬ 
bution and only allow differences in the phase 

factor. The Fourier coefficients can therefore be written 
as the product of a real function /(w) defining the spec¬ 
tral distribution (independent of 7) and a phase factor 
(/)G)(a;) (dependent on 7): 

^^(a;) =/(a;)e*'^‘^’(“). ( 9 ) 

The field auto-correlation function averaged over the en¬ 
semble of realizations then becomes: 


0 0 


X 



( 10 ) 


We choose the phase factor with a linear dependence on 
frequency, and interpret the ensemble averaging as a con¬ 
tinuous sum over the phases, i.e. 


(j)^'^\ui) —>- 70 ; 



(lla) 

(llb) 


so that the last line in (10) yields a delta function and re¬ 
duces the two frequency integrals to a single one. Asking 
that (10) matches the thermal-light first-order correla¬ 
tion function (u) restricts the lineshape /(w) so that the 
electric field for each individual realization (8) becomes: 

= \lnih ( 12 ) 


where the first-order correlation function between the i 
and j Cartesian components of the field, for thermal light, 
is [ 16 ]: 

fc ^00 

with n{io) = _ i)-i being the average pho¬ 

ton number at temperature T and ks the Boltzmann’s 
constant. 

Requiring that ( 4 ) describes the same molecular dy¬ 
namics as that after thermal-light excitation implies a 
matching of first-order correlation functions according to: 

(A(^)(r2)*F;('^)(n))^ = - T 2 ). ( 7 ) 


Using this field as the individual field realizations, the 
individual excitation events (2) become 

Nc 

= K'£paCa{t,j)e-^--*\a), ( 13 ) 

a=l 

where the complex constant K and coefficients Ca{t,'y) 
are detailed in the Supplementary Information (SI). The 
state represents a contribution to the total den¬ 

sity matrix following excitation from one single realiza¬ 
tion of the field. The total density matrix (o) is then 
obtained by averaging over the ensemble according to 
our definition (11b). 

It is easy to analytically show that, with our definition 
of the field realizations (12) and using (11b), the density 
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FIG. 1: Sample of the molecular (a) population paa and (b) 
coherence pap' the total response, calculated from the sum 
over the single excitation events (15) —lines, matches 
the response to thermal-light excitation ( )— p^^{t) dots. 
Note the different scales and amplitudes of coherence vs. pop¬ 
ulation. Parameters used are representative of a photosyn¬ 
thetic light-harvesting complex: £„ = 1 eV, ep = 1.1 eV, 
Pa = Pp = 1 e.A. Convergence of is obtained using 

20 t realizations, where t is the time in femtoseconds. 

matrix (o) built from the ensemble of individual excita¬ 
tion events recovers that obtained directly from thermal 
light excitation (5), i.e. 

= (14) 

Specifically, ( ) becomes 

/ -l-oo 

d7Ca(t,7)C^(i,7)- 

-oo 

(15) 

Sample responses are presented in Fig. (1), showing the 
perfect match and equivalence between methods. As in 
[15, 17, 18], where the incoherent radiation is modeled as 
white-noise, we find that sudden turn-on incoherent exci¬ 
tation can still lead to coherent excitation, with the am¬ 
plitude of the coherence inversely proportional to the en¬ 
ergy splitting between eigenstates. In our approach, the 
coherence survives destructive interference between the 
individual excitation events with primary contribution 
from realizations with 0 < 7 < f. Although non-zero, 
the coherences are of constant amplitude and quickly, i.e. 
after a few femtoseconds, become negligible compared to 
linearly growing populations. 

Interpretations — By construction, the spectral den¬ 
sity of each field realization (12) is identical to that of 
Planck’s radiation law because of the matched correlation 
functions (7). Hence the field realizations’s coherence 
time matches that of thermal light (about 1.3 fs)—since, 
we recall, the coherence time is defined as the inverse of 
the spectral density FWHM and therefore independent 
of the phase factor, cf. (b). As seen in Fig. (2), the lin¬ 
ear dependence of the phase on frequency (lla) localizes 



FIG. 2: Illustration of (a) the phase and (b) the real 
part of the electric field’s positive frequency component 
( 12 ) representing a sample of field realizations—for 7 = 
{0; 50; 100; 150; 200}. All fields have the same spectral dis¬ 
tribution and coherence time, matching that of blackbody 
radiation. The linear phase (lla) results in localized pulses, 
with a duration of ~5 fs, and acts such that the field realiza¬ 
tions only differ from one another by a translation in time. 


the individual realizations, forming pulses. In addition, 
the Xuj phase factor acts such that the field realizations 
only differ from one another by a time translation—also 
evident from the time representation in Fig. (2). Thus, 
the averaging process over all realizations (lib) becomes 
equivalent to averaging over all times, what recovers the 
classical time average. 

Fig. (o) shows the individual excitation events (13) 
corresponding to the particular field realizations sampled 
in Fig. (2). Each event corresponds to a response from 
an impulsive driving force, i.e. with coherent evolution 
and constant population after excitation [19]. Although 
the field realizations are all identical to within a time 
factor, the individual molecular responses are more com¬ 
plex. Specifically, the field realizations differ by a phase 
factor, reflected in a time shift of the individual diago¬ 
nal elements contributing to the total populations (t). 
On the other hand, coherence-wise, the individual molec¬ 
ular contributions p^°p{t) present, in addition to the time 
translation, variations in amplitude. Consequently, al¬ 
though the individual molecular responses are not phase 
related (because of different onsets), for sudden turn-on, 
the total contribution can still exhibit some non-zero os¬ 
cillatory coherences surviving destructive interferences. 
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(a) 



FIG. 3: Illustration of the individual contributions to the to¬ 
tal molecular response, each correspond to excitation from the 
particular field realizations presented in Fig. (z,); (a) popu¬ 
lation \{a\ip^a\t)}\'^ and (b) coherence 


Discussion - The pure states (13) correspond to ex¬ 
citation by the field realizations (12) and it is possi¬ 
ble to build, from these individual responses, a molec¬ 
ular density matrix equivalent to that resulting from 
thermal-light excitation. Experimentally, this would re¬ 
quire recording the molecular response after excitation 
at all times and averaging over the ensemble of these in¬ 
dividual responses. 

It is noteworthy that the ensemble average can be 
performed in many ways. Similar results fulfilling (14) 
can be obtained with individual realizations that dif¬ 
fer from (12) —e.g. non-pulsed light could be used in¬ 
stead. More generally, the only requirement to achieve 
the thermal result is to match the thermal-light first- 
order correlation function, which requires: (i) match¬ 
ing the spectral distribution (and hence the coherence 
time), and (ii) that the ensemble of field realizations ful¬ 
fills (X S{uj-uj'). These conditions 

will reproduce the stochastic character of the first-order 
correlation of thermal-light. The latter requirement re¬ 
sembles the condition that the Fourier components of a 
stationary random process belonging to different frequen¬ 
cies are uncorrelated [1]. While the Fourier components 
of the derived fields (12) are correlated at the level of the 
individual realizations, this correlation is lost upon aver¬ 
aging over an infinity of independent realizations. This 
is in accordance with the fact that, for a stationary and 


ergodic ensemble of random functions, one can replace 
the auto-correlation function defined by the time aver¬ 
age by that defined as an ensemble average [20]. Our 
results demonstrate that this not only holds at the level 
of the field, but also applies to describe linear excitation 
of molecules. 

We emphasize that our results apply to linear inter¬ 
action, relevant to natural incoherent light. In general, 
retrieving the state of thermal light (i.e. its full den¬ 
sity matrix and match all its properties) from an ensem¬ 
ble of single broadband coherent states (representative 
of pulses) is not possible and already fails for the sim¬ 
plest second-order correlation describing broadband de¬ 
layed excitation [10]. Further, we emphasize that the 
individual realizations have no physical meaning. It is 
only the ensemble average that relates to the natural in¬ 
coherent excitation process. 

In conclusion, the implications of our results are two¬ 
fold for the modelling of linear light-matter interaction: 
conceptually, we showed that linear excitation by thermal 
light can be represented by an ensemble average of deter¬ 
ministic pulse realizations, with no random character— 
all effects of incoherence are recast in the post excitation 
averaging. Practically, we provide, with Eq. (13), an 
alternative description of the molecular excitation that 
allows one to study the dynamics of a molecule at the 
wave function level, which presents great simplifications 
over the conventional density matrix formalism. Further, 
the results motivate experimental studies of natural inco¬ 
herent excitation using pulsed laser techniques. As such, 
our findings are expected to facilitate the study of the 
molecular dynamics after thermal-light excitation, such 
as photosynthetic light-harvesting complexes under nat¬ 
ural excitation conditions. 
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Dynamics of a molecular system excited by one realization of the field 


Here we provide details of the wave function dynamics after excitation by one realization of the field (2). Taking, 
as exciting held, the individual realizations derived in (12) yields: 




(16) 


h V l2TT^eoC^ 

Here, we have assumed the held to lie along the orientation of the transition dipole moment Hq,, and we have dehned 


poo 

p^ir) = / duj u)'^^LO n{uj)e~ 

Jo 


Expanding = e '‘“/^''b^( 1 — e '‘“/'“b^) in terms of power of exponentials yields to a series of terms for 

which the integral can be evaluated, i.e. 


Po(t) 

1=0 


2 (-1)^ 


1=0 


doj ' 


(-l)^V^ 


^-l(^/2+j)fl/ksT+ir]ui 


j / 2[(V2 + j)VfcBT + ir]3/2 ■ 


Note that Pvii) = (~l/*)‘^^Po(i)- We therefore take the time derivative to obtain 


S-v/tt 


00/1 
2 




(-1)-^ [(V2 + j)V'=bT + it] 


- 5/2 


We add the integration in time and dehne: 


1=0 ^ ^ 


dr 


=E 


00/1 
2 


(-IF 


0 [(V2 + j)V^BT + l(r - 7)]®^^ 

h, [Ja,j (t, 7) - Ja,] (0, 7)] - - F yj UJahj { 0 , j) 


(17) 


with F{x) = e fg dy is the Dawson function and: 


= (1/2 +j)VfcBr + /(t- 7 ) 


„itu}„ (t, 7) + 1 


Although (17) involves an inhnite number of terms, the sum converges quickly and, at room temperature, results are 
typically within 2% accuracy accounting for the hrst three terms only. 

Taking 


K = -i 


647r2eofic3 ’ 


( 18 ) 


the dynamics of system (16) reduces to the expression given in the main text (18). 
















